Interactive pseudo-3D visualisation with Spatial Analysis#
1. produce an [interactive pseudo-3D Building Model visualization] - via pydeck - which a user can navigate, query, share that;
i) [colour buildings by type] (to easily visualize building stock)
ii) [includes additional features] (parks, bus rapid transit route, etc.)2. allow the user to execute an application of Spatial Data Science
i) [population estimation] –with a previous census metric population growth rate and projected (future) population are also possible and
ii) a measure of [Building Volume per Capita].3. propose several [Geography and Sustainable Development Education conversation starters] for Secondary and Tertiary level students
The village processing option is meant for areas with no more than for 2 500 buildings.
#load the magic
%matplotlib inline
import os
from pathlib import Path
import requests
import overpass
import osm2geojson
import numpy as np
import pandas as pd
import geopandas as gpd
import shapely
from shapely.geometry import Polygon, shape, mapping
#shapely.speedups.disable()
import json
import geojson
import topojson as tp
import fiona
import city3D
import matplotlib.pyplot as plt
import pydeck as pdk
1. Interactive Visualization#
Harvest OpenStreetMap - Query the Overpass API from within Jupyter and convert to .geojson.
This is done: large area -> focus area or State (Province) -> Village (neighourhood / campus)
large = 'Western Cape'
focus = 'Mamre'
osm_type = 'relation'
query = """
[out:json][timeout:30];
// --when areas have duplicate names given the world has a limited amount of uniquely named places
area[name='{0}'] ->.b;
// -- target area ~ can be way or relation
wr(area.b)[name='{1}'];
map_to_area -> .a;
// I want all buildings
(way['building'](area.a);
// and relation type=multipolygon ~ to removed courtyards from buildings
relation["building"]["type"="multipolygon"](area.a);
);
out count;
out geom 2500;
//out body;
//>;
//out skel qt;
""".format(large, focus)#osm_type, focus)
url = "http://overpass-api.de/api/interpreter"
r = requests.get(url, params={'data': query})
#rr = r.read()
gj = osm2geojson.json2geojson(r.json())
village will return a maximum of 2 500 buildings in any focus area.
#- some print statements
print('')
print(focus, 'has', r.json()['elements'][0]['tags']['ways'], 'buildings')
if int(r.json()['elements'][0]['tags']['ways']) < 2500:
print('\n\033[1mAll the buildings\033[0m in', focus, 'have been harvested')
else:
print('\n', int(r.json()['elements'][0]['tags']['ways'])-2500, "buildings have not been harvested.")
#-- try focus = 'Salt River' to see how a small urban suburb will perform or go over to the Suburb folder
Mamre has 2214 buildings
All the buildings in Mamre have been harvested
Please do not burden the OpenStreetMap server with excessive calls for data.
If you need to investigate a larger area (> 2 500 buildings); choose suburb please.
# have a look at a random building
gj['features'][2]
{'type': 'Feature',
'properties': {'type': 'way',
'id': 328118449,
'tags': {'building': 'yes', 'building:levels': '1'},
'nodes': [3349433548, 3349433562, 3349433561, 3349433547, 3349433548]},
'geometry': {'type': 'Polygon',
'coordinates': [[[18.4705985, -33.5053849],
[18.4706523, -33.5056294],
[18.4705655, -33.5056426],
[18.4705116, -33.5053982],
[18.4705985, -33.5053849]]]}}
We assume a building level is 2.8 meters high and add another 1.3 meters (to account for the roof) and create a new attribute height.
The Python code to execute the .calc_Bldheight function is in the city3D.py script
# -- execute function. calculate building height
#city3D.calc_Bldheight(gj)
#city3D.calc_Bldheight(gj, is_geojson=True)
city3D.calc_Bldheight(gj, is_geojson=True, output_file='./data/fp_j.geojson')
pydeck needs an area and a center, harvest from osm.
query = """[out:json][timeout:30];
area[boundary=administrative][name='{0}'] -> .a;
(
way[amenity='university'][name='{2}'](area.a);
{1}[place][place~"sub|town|city|count|state|village|borough|quarter|neighbourhood"][name='{2}'](area.a);
);
out geom;
""".format(large, osm_type, focus)
url = "http://overpass-api.de/api/interpreter"
r = requests.get(url, params={'data': query})
#rr = r.read()
area = osm2geojson.json2geojson(r.json())
#area['features'][0]
#read into .gpd
gdf = gpd.GeoDataFrame.from_features(area['features'])
#-- some relation aoi's are many relations ~ extract the 'place'
if osm_type == 'relation' and len(gdf) > 1:
gdf.dropna(subset = ["tags"], inplace=True)
for i, row in gdf.iterrows():
if row.tags != None and row.tags != np.nan and 'place' in row.tags:
focus = row
trim = pd.DataFrame(focus)
trim = trim.T
gdf = gpd.GeoDataFrame(trim, geometry = trim['geometry'])
#gdf = gdf.set_crs(4326)
# -- get the location for pydeck
[xy] = gdf.geometry.centroid
bbox = [gdf.total_bounds[0], gdf.total_bounds[1],
gdf.total_bounds[2], gdf.total_bounds[3]]
#xy.x
#gdf
~ In order to make the most of the semantic data we need to extract the osm_tags from the dictionary: and add it as tooltips to the visualization.
data = './data/fp_j.geojson'
json = pd.read_json(data)
build_df = pd.DataFrame()
# Parse the geometry out to geopandas
gdf = gpd.GeoDataFrame.from_features(json['features'])
build_df["height"] = round(json["features"].apply(lambda row: row["properties"]["building_height"]), 1)
build_df["plus_codes"] = json["features"].apply(lambda row: row["properties"]["plus_code"])
#we want to display data so extract values from the dictionary
build_df["id"] = json["features"].apply(lambda row: row["properties"]["osm_id"])
build_df["tags"] = json["features"].apply(lambda row: row["properties"])#["osm_tags"])
build_df['levels'] = build_df['tags'].apply(lambda x: x.get('building:levels')).astype(float)
build_df['building'] = build_df['tags'].apply(lambda x: x.get('building'))
build_df['building:use'] = build_df['tags'].apply(lambda x: x.get('building:use'))
build_df['address'] = build_df['tags'].apply(lambda x: x.get('address'))
build_df['building:flats'] = build_df['tags'].apply(lambda x: x.get('building:flats'))
build_df['building:units'] = build_df['tags'].apply(lambda x: x.get('building:units'))
build_df['amenity'] = build_df['tags'].apply(lambda x: x.get('amenity'))
build_df['social_facility'] = build_df['tags'].apply(lambda x: x.get('social_facility'))
build_df['residential'] = build_df['tags'].apply(lambda x: x.get('residential'))
build_df['beds'] = build_df['tags'].apply(lambda x: x.get('beds'))
build_df['rooms'] = build_df['tags'].apply(lambda x: x.get('rooms'))
build_df = gpd.GeoDataFrame(build_df, geometry=gdf.geometry, crs=4326)
#- look
build_df.head(2)
| height | plus_codes | id | tags | levels | building | building:use | address | building:flats | building:units | amenity | social_facility | residential | beds | rooms | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 4.1 | 4FRWFFVC+P79 | 328118446 | {'osm_id': 328118446, 'building': 'yes', 'buil... | 1.0 | yes | None | None | None | None | None | None | None | None | None | POLYGON ((18.47060 -33.50579, 18.47070 -33.505... |
| 1 | 6.9 | 4FRWFFVC+H9Q | 328118447 | {'osm_id': 328118447, 'building': 'church', 'b... | 2.0 | church | None | Kerk Street Mamre 7347 Cape Town | None | None | place_of_worship | None | None | None | None | POLYGON ((18.47092 -33.50592, 18.47093 -33.505... |
# have a look at the building type and amenities available
#df2['bld'].unique()
build_df['building'].unique()
array(['yes', 'church', 'house', 'public', 'civic', 'office', 'retail',
'clinic', 'school', 'cabin', 'garage', 'greenhouse', 'roof',
'kindergarten', 'clubhouse', 'service', 'detached', 'shed'],
dtype=object)
len(build_df)
2187
#- some data wrangling to account for when building:use is different from the original purpose
#- (building=warehouse now loft apartments or =church now office, etc.)
df2 = build_df.copy()
df_res = build_df[build_df['building:use'] == 'residential']
#df_res = df2[df2['building:use'] != None]
df_res = df_res[~df_res['building:use'].isna()]
df2.loc[df_res.index, 'building'] = df_res['building:use']
#build_df.plot()
# colour buildings based on use / amenity
def color(bld):
#- formal house
if bld == 'house' or bld == 'semidetached_house' or bld == 'terrace': #- add maisonette, duplex, etc.
return [255, 255, 204] #-grey
if bld == 'apartments':
return [252, 194, 3] #-orange
#- informal structure / social housing / student
if bld == 'residential' or bld == 'dormitory' or bld == 'cabin':
return [119, 3, 252] #-purple
if bld == 'garage' or bld == 'parking':
return [3, 132, 252] #-blue
if bld == 'retail' or bld == 'supermarket':
return [253, 141, 60]
if bld == 'office' or bld == 'commercial':
return [185, 206, 37]
if bld == 'school' or bld == 'kindergarten' or bld == 'university' or bld == 'college':
return [128, 0, 38]
if bld == 'clinic' or bld == 'doctors' or bld == 'hospital':
return [89, 182, 178]
if bld == 'community_centre' or bld == 'service' or bld == 'post_office' or bld == 'hall' or bld == 'civic' \
or bld == 'townhall' or bld == 'police' or bld == 'library' or bld == 'fire_station' :
return [181, 182, 89]
if bld == 'warehouse' or bld == 'industrial':
return [193, 255, 193]
if bld == 'hotel':
return [139, 117, 0]
if bld == 'church' or bld == 'mosque' or bld == 'synagogue':
return [225, 225, 51]
else:
return [255, 255, 204]
#build_df["fill_color"] = build_df['combine'].apply(lambda x: color(x))
df2["fill_color"] = df2['building'].apply(lambda x: color(x))
df2.head(2)
#build_df.coordinates[0]
| height | plus_codes | id | tags | levels | building | building:use | address | building:flats | building:units | amenity | social_facility | residential | beds | rooms | geometry | fill_color | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 4.1 | 4FRWFFVC+P79 | 328118446 | {'osm_id': 328118446, 'building': 'yes', 'buil... | 1.0 | yes | None | None | None | None | None | None | None | None | None | POLYGON ((18.47060 -33.50579, 18.47070 -33.505... | [255, 255, 204] |
| 1 | 6.9 | 4FRWFFVC+H9Q | 328118447 | {'osm_id': 328118447, 'building': 'church', 'b... | 2.0 | church | None | Kerk Street Mamre 7347 Cape Town | None | None | place_of_worship | None | None | None | None | POLYGON ((18.47092 -33.50592, 18.47093 -33.505... | [225, 225, 51] |
To show the potential and power of 3D City Models we can add additional features to the visualization; namely: bus rapid transit, parks, agricultural land and waterways (streams). We get this from OpenStreetMap as well..
query = """[out:json][timeout:30];
// --main area
area[name='{0}']->.b;
// -- target area ~ can be way or relation
{1}(area.b)[name='{2}'];
map_to_area -> .a;
(
// query
way["sport"](area.a);
way['leisure'="park"](area.a);
);
// print results
out geom;
//out body;
//>;
//out skel qt;
""".format(large, osm_type, focus)
#url = "http://overpass-api.de/api/interpreter"
p = requests.get(url, params={'data': query})
#rr = r.read()
green_spaces = osm2geojson.json2geojson(p.json())
query = """[out:json][timeout:30];
// --main area
area[name='{0}']->.b;
// -- target area ~ can be way or relation
{1}(area.b)[name='{2}'];
map_to_area -> .a;
(
// query
way['waterway'='stream'](area.a);
);
// print results
out body;
>;
out skel qt;""".format(large, osm_type, focus)
#url = "http://overpass-api.de/api/interpreter"
w = requests.get(url, params={'data': query})
#rr = r.read()
water_spaces = osm2geojson.json2geojson(w.json())
query = """[out:json][timeout:30];
// --main area
area[name='{0}']->.b;
// -- target area ~ can be way or relation
{1}(area.b)[name='{2}'];
map_to_area -> .a;
(
// query
way['landuse'='farmland'](area.a);
);
// print results
out body;
>;
out skel qt;
""".format(large, osm_type, focus)
url = "http://overpass-api.de/api/interpreter"
f = requests.get(url, params={'data': query})
#rr = r.read()
p_spaces = osm2geojson.json2geojson(f.json())
# the bus route ~~ note we only choose routes with a 'colour' tag
query = """
[out:json][timeout:30];
area[name='{0}'];
// -- target area ~ can be way or relation
// gather results
(
// query part for: “"bus route"”
relation["type"="route"]["route"="bus"]['operator'="MyCiTi"]['colour'](area);
);
// print results
out body;
>;
out skel qt;
""".format(large)
url = "http://overpass-api.de/api/interpreter"
r = requests.get(url, params={'data': query})
#rr = r.read()
r_lines = osm2geojson.json2geojson(r.json())
# have a look at a random bus route
r_lines['features'][0]['properties']
{'type': 'relation',
'id': 947075,
'tags': {'colour': '#AACDD2',
'from': 'Civic Centre',
'name': 'Bus A01: Civic Centre – Airport [Suspended]',
'network': 'Cape Town IRT',
'operator': 'MyCiTi',
'public_transport:version': '2',
'ref': 'A01',
'route': 'bus',
'to': 'Airport',
'type': 'route'}}
# extract path and assign colour ~~ so the visualization matches the official documentation
Rgdf = gpd.GeoDataFrame.from_features(r_lines['features'])
Rgdf = Rgdf.explode(index_parts=True)
def coords(geom):
return list(geom.coords)
Rgdf['path'] = Rgdf.apply(lambda row: coords(row.geometry), axis=1)
Rgdf['name'] = Rgdf['tags'].apply(lambda x: x.get('name')\
if isinstance(x, dict) else np.nan)
Rgdf['colour'] = Rgdf['tags'].apply(lambda x: x.get('colour', np.nan)\
if isinstance(x, dict) else np.nan)
Rgdf = Rgdf[Rgdf['colour'].notna()]
def hex_to_rgb(h):
h = h.lstrip("#")
#h = h.replace('#', '')
return tuple(int(h[i : i + 2], 16) for i in (0, 2, 4))
Rgdf["colour"] = Rgdf["colour"].apply(hex_to_rgb)
Rgdf.head(3)
| type | id | tags | nodes | geometry | path | name | colour | ||
|---|---|---|---|---|---|---|---|---|---|
| 0 | 0 | relation | 947075 | {'colour': '#AACDD2', 'from': 'Civic Centre', ... | NaN | LINESTRING (18.42874 -33.91977, 18.42884 -33.9... | [(18.4287357, -33.9197677), (18.4288394, -33.9... | Bus A01: Civic Centre – Airport [Suspended] | (170, 205, 210) |
| 1 | 0 | relation | 947076 | {'colour': '#AACDD2', 'from': 'Airport', 'name... | NaN | LINESTRING (18.56274 -33.96201, 18.56253 -33.9... | [(18.5627383, -33.9620148), (18.5625333, -33.9... | Bus A01: Airport – Civic Centre [Suspended] | (170, 205, 210) |
| 1 | relation | 947076 | {'colour': '#AACDD2', 'from': 'Airport', 'name... | NaN | LINESTRING (18.56357 -33.96227, 18.56274 -33.9... | [(18.563566, -33.9622656), (18.5627383, -33.96... | Bus A01: Airport – Civic Centre [Suspended] | (170, 205, 210) |
## ~ (x, y) - bl, tl, tr, br ~~ or ~~ sw, nw, ne, se
#area = [[[18.4377, -33.9307], [18.4377, -33.9283], [18.4418, -33.9283], [18.4418, -33.9307]]]
area = [[[bbox[0], bbox[1]], [bbox[0], bbox[3]],
[bbox[2], bbox[3]], [bbox[2], bbox[1]]]]
## ~ (y, x)
view_state = pdk.ViewState(latitude=xy.y, longitude=xy.x, zoom=16.5, max_zoom=19, pitch=72,
bearing=80)
land = pdk.Layer(
"PolygonLayer",
#"GeoJsonLayer",
area,
#build_df,
stroked=False,
# processes the data as a flat longitude-latitude pair
get_polygon="-",
get_fill_color=[0, 0, 0, 1],
#material = True,
#shadowEnabled = True
)
building_layer = pdk.Layer(
"PolygonLayer",
#"GeoJsonLayer",
df2,
#id="geojson",
opacity=0.3,
stroked=False,
get_polygon="geometry.coordinates",
filled=True,
extruded=True,
wireframe=False,
get_elevation="height",
#get_fill_color="[255, 255, 255]", #255, 255, 255
get_fill_color="fill_color",
get_line_color="fill_color",#[255, 255, 255],
#material = True,
#shadowEnabled = True,
auto_highlight=True,
pickable=True,
)
greenspaces_layer = pdk.Layer(
"GeoJsonLayer",
green_spaces,
opacity=0.5,
stroked=False,
filled=True,
wireframe=True,
get_fill_color="[14, 140, 58]",
get_line_color='[14, 140, 58]',
)
water_layer = pdk.Layer(
"GeoJsonLayer",
water_spaces,
opacity=0.8,
stroked=False,
filled=True,
wireframe=True,
get_fill_color="[35, 35, 142]",
get_line_color='[35, 35, 142]',
)
r_layer = pdk.Layer(
type="PathLayer",
data=Rgdf,
get_color='colour', #'[245, 51, 58]',
#width_scale=20,
#width_min_pixels=8,
get_path="path",
get_width=5,
auto_highlight=False, # change to True if route query
pickable=False, # change to True if route query
)
p_layer = pdk.Layer(
"GeoJsonLayer",
p_spaces,
opacity=0.2,
stroked=False,
filled=True,
wireframe=True,
get_fill_color="[170, 83, 3]",
get_line_color='[170, 83, 3]',
)
tooltip = {"html": "<b>Levels:</b> {levels} <br/> <b>Address:</b> {address}\
<br/> <b>Plus Code:</b> {plus_codes} <br/> <b>Building Type:</b> {building}"}
#change the tooltip to show bus routes and comment out the previous
#tooltip = {"html": "<b>Route:</b> {name} <br/>"}
r = pdk.Deck(layers=[greenspaces_layer, p_layer, water_layer, r_layer, building_layer],
#views=[{"@@type": "MapView", "controller": True}],
initial_view_state=view_state,
map_style = 'dark_no_labels', #pdk.map_styles.LIGHT,
tooltip=tooltip)
#save
r.to_html("./result/interactiveOnly.html")
#r.show()